#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

// 定义右端项函数
vector<double> f(double t, vector<double>& u, double mu) {
    vector<double> res(6);  // 初始化结果向量
    res[0] = u[3];
    res[1] = u[4];
    res[2] = u[5];
    res[3] = 2*u[4] + u[0] - mu*(u[0]+mu-1)/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu-1, 2), 1.5) - (1-mu)*(u[0]+mu)/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu, 2), 1.5);
    res[4] = 2*u[3] + u[1] - mu*u[1]/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu-1, 2), 1.5) - (1-mu)*u[1]/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu, 2), 1.5);
    res[5] = -mu*u[2]/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu-1, 2), 1.5) - (1-mu)*u[2]/pow(u[1]*u[1]+u[2]*u[2]+pow(u[0]+mu, 2), 1.5);
    return res;
}

// 实现Adams-Bashforth methods
vector<double> AdamsBashforth(double tn, double h, vector<double>& un, double mu, int p) {
    vector<double> u_next(6);  // 初始化下一步解
    if (p == 1) {
        vector<double> f_n = f(tn, un, mu);  // 计算f_n
        for (int i = 0; i < 6; ++i) {
            u_next[i] = un[i] + h*f_n[i];  // 计算下一步解
        }
    } else if (p == 2) {
        vector<double> f_n = f(tn, un, mu);  // 计算f_n
        vector<double> f_n_1 = f(tn-h, un, mu);  // 计算f_{n-1}
        for (int i = 0; i < 6; ++i) {
            u_next[i] = un[i] + h*(1.5*f_n[i] - 0.5*f_n_1[i]);  // 计算下一步解
        }
    } else if (p == 3) {
        vector<double> f_n = f(tn, un, mu);  // 计算f_n
        vector<double> f_n_1 = f(tn-h, un, mu);  // 计算f_{n-1}
        vector<double> f_n_2 = f(tn-2*h, un, mu);  // 计算f_{n-2}
        for (int i = 0; i < 6; ++i) {
            u_next[i] = un[i] + h*(23.0/12*f_n[i] - 4.0/3*f_n_1[i] + 5.0/12*f_n_2[i]);  // 计算下一步解
        }
    } else if (p == 4) {
        vector<double> f_n = f(tn, un, mu);  // 计算f_n
        vector<double> f_n_1 = f(tn-h, un, mu);  // 计算f_{n-1}
        vector<double> f_n_2 = f(tn-2*h, un, mu);  // 计算f_{n-2}
        vector<double> f_n_3 = f(tn-3*h, un, mu);  // 计算f_{n-3}
        for (int i = 0; i < 6; ++i) {
            u_next[i] = un[i] + h*(55.0/24*f_n[i] - 59.0/24*f_n_1[i] + 37.0/24*f_n_2[i] - 3.0/8*f_n_3[i]);  // 计算下一步解
        }
    }
    return u_next;
}

int main() {
    double mu = 0.012277471;  // 常数µ
    double T = 17.065216560157;  // 最终时间
    double h = 0.01;  // 步长
    vector<double> u = {0.994, 0, 0, 0, -2.00158510637908, 0};  // 初始值
    int N = (int)(T/h);  // 时间步数
    vector<vector<double>> sol(N+1, vector<double>(6));  // 初始化解向量
    sol[0] = u;  // 存储初始值
    for (int n = 0; n < N; ++n) {
        double tn = n*h;  // 当前时间
        vector<double> un = sol[n];  // 当前解
        vector<double> u_next = AdamsBashforth(tn, h, un, mu, 4);  // 计算下一步解
        sol[n+1] = u_next;  // 存储下一步解
    }
    // 输出最终时间下的解向量
    for (int i = 0; i < 6; ++i) {
        cout << "u_" << i+1 << "(" << T << ") = " << sol[N][i] << endl;
    }
    return 0;
}